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Abstract 



In this paper we present NAN0TCAD2D, a code for the simulation of 
the electrical properties of semiconductor-based nanoelectronic devices and 
structures in two-dimensional domains. Such code is based on the solution 
of the Poisson/Schrodinger equation with density functional theory and of 
the continuity equation of the ballistic current. NAN0TCAD2D can be ap- 
plied to structures fabricated on III- IV, strained-silicon and silicon-germanium 
heterostructures, CMOS structures, and can easily be extended to new ma- 
terials. In particular, in the case of SiGe heterostructures, it computes the 
effects of strain on the energy band profiles. The effects of interface states at 
the air /semiconductor interfaces, particularly significant in the case of devices 
obtained by selective etching, are also properly taken into account. 

I. PROGRAM SUMMARY 



Title of program: NAN0TCAD2D 

Program summary [//ZL.-:http://www. phantomshub.com 



1 



Program available at: URL:http://www.phaiitomshub.com. (can be freely used via a web 
interface) 

Computer on which the program has been tested: Linux PCs. 
Programming Language used: FORTRAN 77. 
Memory used to execute with typical data: 50-100 Mb. 

Keywords: density functional theory, quantum simulation, ballistic transport, nanoscale 

devices, nanoelectronics, nanotechnology computer aided design. 

Field of application: simulation of nanoscale semiconductor devices and structures. 

Method of solution: Self- consistent solution of the Poisson-Schrodinger equation with density 

functional theory, in the approximation of local density, and of the continuity equation for 

the ballistic current with the Newton-Raphson algorithm. 

Restrictions: The present simulation tool works only with rectangular grids. 



II. INTRODUCTION 

The electrical properties of nanoscale semiconductor devices and structures are typically 
determined by quantum confinement, that strongly affects the density of states of electrons 
and holes, and, by ballistic transport, that takes place if device length is smaller than the 
scattering length. 

An efficient simulation tool, aimed at understanding device behavior and at designing 
optimized structures, must therefore include such aspects, and achieve at the same time a 
reasonable trade-off between efficiency and accuracy. 

Here we present NAN0TCAD2D, a program for the simulation of nanoelectronic semi- 
conductor devices, based on the self-consistent solution of Poisson/Schrodinger equation and 
of the continuity equation for electrons and holes in the case of ballistic transport. The pro- 
gram allows to consider quantum confinement in one or two dimensions, and to address more 
complex structures that can be divided into regions in which different types of confinement 
are present. 
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NAN0TCAD2D is oriented at two main classes of devices: i) quantum wires, i.e., struc- 
tures with translation symmetry in one direction, that are completely described by the 
geometry of the cross section; and ii) ballistic field effect devices, again with translation 
symmetry in a direction perpendicular to that of electron motion. 

Ballistic transport is simply included by assuming that the occupation factor of states in- 
jected from a given reservoir is that of the reservoir itself, i.e., obeys Fermi-Dirac distribution 
with the reservoir's Fermi energy. 

In addition, the code implements a simplified model for localized states at the semi- 
conductor surface exposed to air. Indeed, the depletion of a one-Dimensional or a two- 
Dimensional Electron Gas (IDEG or 2DEG) due to acceptor-like surface states is of primary 
importance in determining the electrical properties of narrow quantum wires or structures 
with large exposed surfaces. 

At present, the code allows to consider common semiconductor materials, such as silicon, 
AlGaAs, InGaAs, strained silicon and silicon germanium. In the case of silicon germanium, 
in particular, we have developed a procedure for taking into account the effect of strain 
caused by different lattice constants on the band structure. 

The paper has the following structure: in section III we describe the physical model 
implemented in the code; in section IV we discuss the numerical aspects algorithms. In 
section V simulation examples are presented and in section VI the conclusions. The appendix 
contains the user's manual. 

III. PHYSICAL MODEL 

The density of states for electron and holes depends on the degree of quantum con- 
finement assumed in the different device regions. There are three possibilities: If quantum 
confinement is strong in both directions, the density of states is obtained by solving the two- 
dimensional Schrodinger equation. If quantum confinement is strong only in one direction 
(say, x) the density of states is written as a sum of two-dimensional subbands, the edges 
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of which are obtained by solving the Schrodinger equation in the x direction for each mesh 
point along the y axis. If confinement can be considered very weak, the density of states of 
the bulk material is used. 

For materials with degenerate or quasi-degenerate minima in the conduction band, such 
as for example silicon, that has six degenerated minima, the density of states is com- 
puted with the effective mass approximation for each minimum, taking into account mass 
anisotropy. The same procedure is applied to degeneracy of valence band maxima. 

In the following, we will discuss in some detail the expressions for the density of states 
and carrier density considering only one band valley, for simplicity of notation. Extension 
to multiple valleys is straightforward. 

A. Charge density for two-dimensional quantum confinement (IDEG-IDHG) 

Let us consider a region where the confinement for electrons is strong in both x and y 
directions. In such a case, the local density of states per unit of volume and energy near a 
conduction band minimum is given by: 



where u{E — E-i) is the Heavyside function, ^'j is the solution of the Schrodinger equation 
in two dimensions, i.e.. 



Ei is the corresponding eigenvalue, Ec is the conduction band, m^, s = x^y^ z is the effective 
mass in the direction denoted by the pedix, fi is the reduced Planck's constant. At this 
point, by integrating the density of states multiplied by the Fermi-Dirac occupation factor, 
the quantum electron density can be expressed as: 



NiD{E,x,y) 




(1) 






(3) 
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where -F_i/2 is the Fermi-Dirac integral of order —1/2 and Ef is the Fermi energy [5]. 

In order to compute the hole concentration we have to solve Schrodinger equation for 
heavy holes and for light holes. Therefore, the conduction band in (2) is substituted by the 
inverted valence band — and the eigenvalues —E^, are obtained. 

Therefore the hole concentration p becomes: 

where is the effective mass for holes in the z direction. 

B. One-dimensional quantum confinement: two-dimensional electron or hole gas 

(2DEG-2DHG) 

In the case of strong confinement in only one direction (for example along the x-direction) , 
we assume that the density of states can be decomposed in a quantum term along the 
confined direction [x) and a semiclassical term in the other directions. The one-dimensional 
Schrodinger equation for electrons in the x direction for a mesh point y can be written as: 



^ :<^>i + E,{x,y)^i^Ei{y)^i (5) 



2mx dx^ 

As a consequence, the available states for electrons are grouped into two-dimensional 
subbands and the density of states can be expressed as follows: 



/tTL 777, 

N2d{E, X, y) - Yl y) I' - ^^(^)) (6) 



2TTh 

The electron density therefore is: 

kBT^myfri; 



n 



Similar considerations apply to holes. 



1 -I- exp 



Ef — Ej 



(7) 
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C. Effects of surface states 



In the simulation of narrow semiconductor devices obtained with selective etching, the 
effects of states at the exposed surface are very important and must be taken into account 
in order to reproduce with accuracy the experimental results. In particular, these states can 
act as donors or acceptors and hence deeply affect the carrier distribution within the device. 

In order to correctly model the phenomenon, we have used a simple model based on 
two parameters that is typically applied to metal-semiconductor contacts [6] and has been 
recently validated for air-semiconductor interfaces [7]. 

The two parameters are the density of interface states per unit energy per unit area Ds 
[eV~^ cm~^] and the energy difference between the vacuum level Eo and the Fermi energy 
that ensures a neutral charge at the interface. States with energy below Eg — are donors 
and states with higher energy are acceptors. 

Surface charge per unit surface can then be expressed as Qg = —qDg [Ef — {Eg — $*)], 
where — g is the electron charge. 



All charge concentrations considered represent the source term of the two-dimensional 
Poisson equation: 



where e is the dielectric constant, q is the electron charge, Ps is the term of surface charge per 
unit volume, and the ionized donor and acceptor concentrations, respectively [6]. 
Energy bands depend on the potential as: 



D. Poisson Equation 



V • (eV$) 



q [-n[$] + + iV+[$] - Ar-[$] + 



(8) 



Ec{x,y) = Ec{x,y)\^=o - q^x,y) 



Ev{x,y) = Ey{x,y)\^=o - q^{x,y). 



(9) 
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Potential and charge density profiles in equilibrium are computed by solving the set of 
non linear partial differential equations described above. The case of ballistic transport is 
examined in the next section. 

E. Ballistic transport for electrons and holes. 

When carriers are injected into a semiconductor device, they are hkely to be scattered 
by a number of possible sources, including acoustic and optical phonons, ionized impurities, 
defects, interfaces and other carriers. If, however, device length is smaller than the mean free 
path it is very likely for carriers to traverse the device without suffering scattering events. 
Our code includes this type of "ballistic" transport. 

Ballistic transport is implemented here only in the case of one-dimensional quantum 
confinement, when a description in terms of two dimensional subbands is used. Indeed, we 
have shown that such assumption involves a negligible error also in devices with channel 
length of 25 nm [8]. 

Let us consider Figure 1, representing a subband profile for electrons along the channel. 
Carriers are injected into the channel from a reservoir (source) and contribute to the current 
only if they overcome the barrier modulated by the gate voltage and, to a lesser degree, by 
the drain voltage (DIBL). 

If we neglect the interaction between electrons and ions and among electrons, we can 
simply assume that electrons with injected longitudinal energy lower than the subband 
maximum are refiected back to their originating contact, while the others are transmitted 
over the barrier and contribute to the current [9,10]. 

Therefore, for each subband we evaluate the subband maximum Eimax and the corre- 
sponding longitudinal position Umax- AH electrons with longitudinal energy lower than Eimax 
are in equilibrium with the originating contact, while electrons with longitudinal energy 
higher than Eimax conserve the chemical potential of the injecting reservoir. The occupation 
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factor / is therefore: 



f{E,EF)={ 



1 + 
1 + 



l + e^P(^ 



if y > ?/max, E < Ei, 

if E > E- 



(10) 



where -EiT^^ {Efd) is the source (drain) Fermi energy. If we write the total energy E as 



h k '~ 

E — Ey-\- Ez, where the term Ey — Ei + is the longitudinal energy, and E^ — fi^k1/2m 



21,2, 



the transverse energy, the density of states reads: 



N^D{Ey,Ez) dE = 2 ^1^, 



2 ^niym^ 

h? yjEyEz 



dEy dEz 



(11) 



and the electron density is accordingly given by: 
Jo 



^ ' Ey--2dEy / E-y{E,EF)dEz 



h? 



(12) 



The electron current is evaluated assuming that there is no tunnel current through the 
barrier so that only the electrons with longitudinal energy higher than Eimax can contribute. 

^ y^myUiz 1 



Jn 



X 



dE, 
1 



1 



2E„ 



_l + exp 1 + exp _ 



dE„ 



(13) 



Similar considerations apply to holes. 



IV. NUMERICAL ASPECTS 

The flow diagram of the algorithm implemented is shown in Figure 2. The program, 
using an initial guess for the potential, starts with a semiclassical solution of the Poisson 
equation. The equation is discretized with the Box-integration method and solved with 
the Newton- Raphson algorithm. Dirichlet boundary conditions are enforced on each metal 
gate and homogeneous Neumann conditions on the rest of the domain boundary. The 
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solution obtained is used as an initial guess for the quantum calculation, where the nested 
Poisson/Schrodinger equation must be solved. 

In order not to degrade convergence speed of the algorithm when also the Schrodinger 
equation has to be solved, we have implemented a simplified version of the predictor-corrector 
scheme proposed in Ref . [1 1] . In this way, instead of solving both equations at each Newton- 
Raphson step, we evaluate eigenfunctions and eigenvalues only at the beginning of a Newton- 
Raphson cycle: for the whole cycle eigenfunctions and the difference between the eigenvalues 
and the energy bands in each point of the domain are assumed to be constant. When a 
Newton-Raphson cycle ends the Schrodinger equation is solved again and a new cycle is 
started. The program ends when the difference between the two-norm of the potential at 
the end of two successive Newton-Raphson cycles is lower then a fixed tolerance. [12] 

In Figure 3(left) a rectangular uniform mesh is shown, where it is possible to notice to 
the subdomain Dij (dashed line) associated to the generic grid point In our code, 

we have distinguished the properties of the structure in point characteristics and material 
characteristics. 

Point characteristics are, for example, the potential, the Fermi level, charge concentra- 
tions, while material properties (e.g., energy gap, electron affinity, etc.) belong to the second 
group. In the figure, is also shown (solid line) the material element connected to the generic 
(i, j) grid point. 

In each subdomain, Poisson and Schrodinger equations are discretized with box integra- 
tion [13] 



A. Poisson Equation 

In particular, for the Poisson equation, after integrating both members over the domain 
Dij, we obtain: 

/ / V ■{eV^x,y))dxdy = - / p{x,y)dxdy, (14) 
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which is discretized as 
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+ 



[-g (p., - n,, + - Ar-^,J + p,,,] ^ "^^^^ ^ ^"^^ (15) 

where /i, = Xj+i — Xj, /cj = y^+i — y^, and pedices i, j denote the quantity in position {xi, yj). 

Proper boundary conditions must be enforced to the equation, such as Dirichlet boundary 
conditions on each metal gate and homogeneous Neumann conditions on the rest of the 
domain boundary. In the first case, the potential $ is fixed, while in the second case, the 
electric field V$ • n = — £ is fixed. 

In the simulation code, we have chosen as reference level for energies the vacuum level 
Eo and hence the potential for each gate point is obtained as: 

$(i,j)=£^o-Co.fe-^F (16) 

where Ep and (f)^^^ represent the Fermi level and the work function of the n-th gate, 
respectively. 

Let us consider a point (i,j) on the boundary and let us make reference to the Fig- 
ure 3(right): observe how the region Dij is much smaller with respect to the case of internal 
point {Dij becomes a quarter in the case of each of four vertexes grid points). In this point, 
we enforce Neumann condition and the discretization of the Poisson equation becomes: 
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[-Q {Pi 



n. 



+ 



k 



'J 



hi-i 



hi-1 



(17) 



B. Schrodinger Equation 



The two-dimensional single-particle Schrodinger equation for electrons, given a conduc- 
tion band profile Ec{x,y), reads: 



(18) 



where ^„(a;, y) represents the n-th eigenfunction, En is the n-th eigenenergy, m is the electron 
effective mass tensor in the plane perpendicular to the direction of propagation, 



m 



rrix 



m,, 



(19) 



In our simulations, we have discarded the exchange- correlation term, since it provides a very 
small contribution. Dirichlet boundary conditions are enforced on the quantum simulation 
domain. [15] 

With the box integration method, we obtain: 
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m 



{hi + hi_i){kj + kj_i] 



yi,j yi—i,j 
{hi + hi_i){kj + %_i) 



(20) 



C. Numerical routines and performance 

The discretized non linear Poisson equation is solved with the Newton-Raphson 
(NR) algorithm. The sparse system of linear algebraic equations of each NR step is solved 
with the package Y12MAF [14], which is based on Gaussian ehmination 

The eigenvalue problem resulting from the discretization of the Schrodinger equation in 
one dimension is solved with the routine TQLI [14], while in two dimensions is solved with 
the method proposed in Ref. [15] that allows to reduce computing time without significant 
losses in accuracy, by solving the problem in the momentum space. The method can be 
applied to structures with inhomogeneous effective mass and can easily be extended to the 
full band structure. 

The computing time on an 1800 MHz Pentium IV CPU strongly depends on the type of 
simulation: In the case of quantum confinement in one direction the CPU running time for 
a 128x61 grid is 37.14 sec, while in the case of quantum confinement in two directions, with 
a 113x148 point grid, the CPU running time is about 95.44 sec. These results represent 
the worst case, with no initial guess of the unknown potential. Finally, in the case of a 
simulation of ballistic current the running time is strongly affected by the initial guess of the 
potential and is between a few minutes and an hour. The initial guess is also very important 
in order to avoid convergence problems of the algorithm. 
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V. EXAMPLES OF SIMULATION 



In this section, two examples of simulations computed on nanoscale devices are shown. 
The first structure is a silicon-germanium high mobility electron waveguide schematically 
represented in Figure 4. It consists of a Sio.8Geo.2 virtual substrate, an 11 nm strained-silicon 
layer in which the IDEG forms, a 5.7 nm undoped Sio.8Geo.2 spacer layer, a 5.7 nm Sio.8Geo.2 
doped layer, with A^^ = 10^^ cm a 35 nm undoped Sio.8Gco.2 spacer and a 15 nm undoped 
silicon cap layer. The second spacer is rather thick, in order to prevent the formation of 
another electron channel in the silicon cap layer. The waveguide is 160 nm wide. Finally, we 
assume that a triple metal gate is deposited over the structure forming a Schottky contact. 
For the purpose of our simulation, the Schottky junction is reverse-biased and assumed to 
be perfectly insulating. 

Quantum confinement of carriers in the horizontal (y) direction is provided by selective 
etching and by the depletion region induced by acceptor states at the exposed surfaces. 
This last effect causes the electrical width of wire to be significantly smaller then the etched 
width. Along the growth {x) direction, as a consequence of the band alignment between 
strained-silicon and silicon-germanium, a quantum well for electrons forms. In particular, 
the strained silicon channel is grown under a tensile strain and thus two valleys of the 
conduction band, along k^., are lowered in energy while the other four valleys are raised. 
This condition is required in order to obtain a confinement region for electrons. In addition, 
only the two lowest conduction band valleys are occupied and the energy splitting between 
valleys ( 120 meV) leads to strongly suppressed intervalley scattering. 

The self-consistent Poisson/Schrodinger equation is discretized onto a nonuniform rect- 
angular grid of 108 X 137 points and the electron concentration has been calculated by solving 
the Schrodinger equation inside the strained silicon channel. Quantum electron density is 
represented in Figure 5, where it is possible to observe how the electrical waveguide width 
is about 95 nm, instead of 160 nm, because of the electron depletion induced by interface 
states at the exposed surfaces. By tuning the voltage applied to the gate, it is possible to 
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vary the electron density in the channel and hence the number of occupied states. Thus, 
referring to the Landauer formula of the quantized conductance G = N (2e^/ h) with N equal 
to the number of propagating modes, it is possible to vary, as a function of applied voltage, 
the quantum conductance in the channel, as shown in Figure 6. 

The second simulation example refers to a "well tempered" ballistic MOSFET with 
channel length of 25nm, proposed by Antoniadis et al. [16] and schematically represented in 
the inset of Figure 7. The oxide thickness is 1.5 nm and the polysihcon gate has a donor 
concentration of 5 x 10^° cm~^. In order to reduce short channel effects a super-halo doping 
is implanted in the channel. The analytic doping profile can be found in Ref. [16]. 

Assuming fully ballistic transport within the channel, we have computed the source- 
to-drain current. The energy barrier that electrons encounter traveling from the source 
towards the drain is represented in Figure 7 for a gate voltage Vqs = 1 V and drain-to- 
source voltage Vds = V. Only the first five subbands are shown. Quantum tunneling has 
not be considered in our model and therefore the transmission coefficient is unity above the 
peak of the barrier and zero below. Thus, only electrons with energy higher than the peak 
can traverse the channel without energy loss and contribute to the total current. 

The transfer characteristics obtained with the MEDICI simulator and with NANOT- 
CAD2D are compared in Figure 8. 

The two-dimensional simulator is available, after registration, at the URL: 
http://www.phantomshub.com. In the directory NAN0TCAD2D it is possible to find 
additional examples of nanoscale semiconductor devices, among which users can find two 
different type high mobility electron waveguides defined by selective etching on a SiGe het- 
erostructures, a nanoscale baUistic sihcon MOSFET and a nanoscale ballistic AlGaAs field 
effect transistor (Figure 9). The parameters of simulation can be varied according to the 
specific user requirements. In particular, input files can be modified by user directly via a 
web interface as an HTML form or uploaded from an external source (Figure 10). In the 
same directory it is also possible to find a detailed tutorial regarding the input data files 
structure. 
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Finally, the hub provides basic visualization capabihties for 2D and 3D plots controlled 
via a web interface (Figure 11). 

VI. DISCUSSION AND FUTURE DEVELOPMENTS 

In this paper we have presented a two dimensional quantum simulator based on the 
solution of the Poisson/Schrodinger equation with the Box- Integration method. The code 
solves also the continuity equation for electrons and holes in the case of ballistic transport, 
where propagating states are populated according to the occupation factor of the originating 
reservoir. NAN0TCAD2D allows to simulate most common semiconductors, such as silicon, 
AlGaAs, InGaAs, strained silicon and silicon germanium and has been successfully used to 
simulate III-IY heterostructures, strained-silicon and silicon germanium heterostructures, 
CMOS structures. In the case of silicon-germanium based devices a dedicated procedure 
allows to take into account the effects of strain on the energy bands and on band alignment. 

The code is presently being improved with full two-dimensional quantum transport and 
quantum tunneling, that becomes relevant for devices with channel length close to 10 nm. 
In addition, we are developing a model for quasi-ballistic transport, which accounts for the 
possibility that a fraction of carriers undergo elastic scattering, that would allow a more 
accurate simulation of nanoscale field effect transistor at room temperature. 
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FIGURES 
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FIG. 1. Energy profile along the channel for a nanoscale FET. Only electrons with energy 
higher than the barrier peak contribute to the current. 
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FIG. 2. Flow diagram of the algorithm implemented. 
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FIG. 3. Dij represents the region associated to the grid point In figure is presented the 

case of an internal point (left) and a point in the edges (right). 




FIG. 4. Silicon-germanium electron waveguides 




FIG. 5. Quantum electron density in the strained silicon channel. 
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FIG. 6. The voltage applied on the external gate allows us to select the number of propagating 
modes in the waveguide. Hence it's possible to vary the quantized conductance in the Si channel 
G = {le^ /h)N as a function of gate voltage. N is the number of propagating modes. For purpose 
of presentation each curve is shifted by one conductance quantum. 
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FIG. 7. Subband energies for a voltage applied to the gate of IV. Transmission coefficient 
through the barrier is assumed to be unity for electrons with energy higher than the peak of the 
barrier and zero for electron with energy lower than the peak. 
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FIG. 8. Transfercharacteristics of the 25 nm MOSFET computed with MEDICI (left, from Ref. 
Antoniadis) and with NAN0TCAD2D (right). 
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FIG. 9. Each user can use all the programs available on the HUB 
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FIG. 10. The parameters of simulation can be varied in accordance with the specific user 



requirements. 




FIG. 11. The results of simulations can be graphically represented 
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